home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / blas / dgemv.f < prev    next >
Text File  |  1997-01-29  |  7KB  |  262 lines
  1.       SUBROUTINE DGEMV ( TRANS, M, N, ALPHA, A, LDA, X, INCX,
  2.      $                   BETA, Y, INCY )
  3. *     .. Scalar Arguments ..
  4.       DOUBLE PRECISION   ALPHA, BETA
  5.       INTEGER            INCX, INCY, LDA, M, N
  6.       CHARACTER*1        TRANS
  7. *     .. Array Arguments ..
  8.       DOUBLE PRECISION   A( LDA, * ), X( * ), Y( * )
  9. *     ..
  10. *
  11. *  Purpose
  12. *  =======
  13. *
  14. *  DGEMV  performs one of the matrix-vector operations
  15. *
  16. *     y := alpha*A*x + beta*y,   or   y := alpha*A'*x + beta*y,
  17. *
  18. *  where alpha and beta are scalars, x and y are vectors and A is an
  19. *  m by n matrix.
  20. *
  21. *  Parameters
  22. *  ==========
  23. *
  24. *  TRANS  - CHARACTER*1.
  25. *           On entry, TRANS specifies the operation to be performed as
  26. *           follows:
  27. *
  28. *              TRANS = 'N' or 'n'   y := alpha*A*x + beta*y.
  29. *
  30. *              TRANS = 'T' or 't'   y := alpha*A'*x + beta*y.
  31. *
  32. *              TRANS = 'C' or 'c'   y := alpha*A'*x + beta*y.
  33. *
  34. *           Unchanged on exit.
  35. *
  36. *  M      - INTEGER.
  37. *           On entry, M specifies the number of rows of the matrix A.
  38. *           M must be at least zero.
  39. *           Unchanged on exit.
  40. *
  41. *  N      - INTEGER.
  42. *           On entry, N specifies the number of columns of the matrix A.
  43. *           N must be at least zero.
  44. *           Unchanged on exit.
  45. *
  46. *  ALPHA  - DOUBLE PRECISION.
  47. *           On entry, ALPHA specifies the scalar alpha.
  48. *           Unchanged on exit.
  49. *
  50. *  A      - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
  51. *           Before entry, the leading m by n part of the array A must
  52. *           contain the matrix of coefficients.
  53. *           Unchanged on exit.
  54. *
  55. *  LDA    - INTEGER.
  56. *           On entry, LDA specifies the first dimension of A as declared
  57. *           in the calling (sub) program. LDA must be at least
  58. *           max( 1, m ).
  59. *           Unchanged on exit.
  60. *
  61. *  X      - DOUBLE PRECISION array of DIMENSION at least
  62. *           ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n'
  63. *           and at least
  64. *           ( 1 + ( m - 1 )*abs( INCX ) ) otherwise.
  65. *           Before entry, the incremented array X must contain the
  66. *           vector x.
  67. *           Unchanged on exit.
  68. *
  69. *  INCX   - INTEGER.
  70. *           On entry, INCX specifies the increment for the elements of
  71. *           X. INCX must not be zero.
  72. *           Unchanged on exit.
  73. *
  74. *  BETA   - DOUBLE PRECISION.
  75. *           On entry, BETA specifies the scalar beta. When BETA is
  76. *           supplied as zero then Y need not be set on input.
  77. *           Unchanged on exit.
  78. *
  79. *  Y      - DOUBLE PRECISION array of DIMENSION at least
  80. *           ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n'
  81. *           and at least
  82. *           ( 1 + ( n - 1 )*abs( INCY ) ) otherwise.
  83. *           Before entry with BETA non-zero, the incremented array Y
  84. *           must contain the vector y. On exit, Y is overwritten by the
  85. *           updated vector y.
  86. *
  87. *  INCY   - INTEGER.
  88. *           On entry, INCY specifies the increment for the elements of
  89. *           Y. INCY must not be zero.
  90. *           Unchanged on exit.
  91. *
  92. *
  93. *  Level 2 Blas routine.
  94. *
  95. *  -- Written on 22-October-1986.
  96. *     Jack Dongarra, Argonne National Lab.
  97. *     Jeremy Du Croz, Nag Central Office.
  98. *     Sven Hammarling, Nag Central Office.
  99. *     Richard Hanson, Sandia National Labs.
  100. *
  101. *
  102. *     .. Parameters ..
  103.       DOUBLE PRECISION   ONE         , ZERO
  104.       PARAMETER        ( ONE = 1.0D+0, ZERO = 0.0D+0 )
  105. *     .. Local Scalars ..
  106.       DOUBLE PRECISION   TEMP
  107.       INTEGER            I, INFO, IX, IY, J, JX, JY, KX, KY, LENX, LENY
  108. *     .. External Functions ..
  109.       LOGICAL            LSAME
  110.       EXTERNAL           LSAME
  111. *     .. External Subroutines ..
  112.       EXTERNAL           XERBLA
  113. *     .. Intrinsic Functions ..
  114.       INTRINSIC          MAX
  115. *     ..
  116. *     .. Executable Statements ..
  117. *
  118. *     Test the input parameters.
  119. *
  120.       INFO = 0
  121.       IF     ( .NOT.LSAME( TRANS, 'N' ).AND.
  122.      $         .NOT.LSAME( TRANS, 'T' ).AND.
  123.      $         .NOT.LSAME( TRANS, 'C' )      )THEN
  124.          INFO = 1
  125.       ELSE IF( M.LT.0 )THEN
  126.          INFO = 2
  127.       ELSE IF( N.LT.0 )THEN
  128.          INFO = 3
  129.       ELSE IF( LDA.LT.MAX( 1, M ) )THEN
  130.          INFO = 6
  131.       ELSE IF( INCX.EQ.0 )THEN
  132.          INFO = 8
  133.       ELSE IF( INCY.EQ.0 )THEN
  134.          INFO = 11
  135.       END IF
  136.       IF( INFO.NE.0 )THEN
  137.          CALL XERBLA( 'DGEMV ', INFO )
  138.          RETURN
  139.       END IF
  140. *
  141. *     Quick return if possible.
  142. *
  143.       IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.
  144.      $    ( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) )
  145.      $   RETURN
  146. *
  147. *     Set  LENX  and  LENY, the lengths of the vectors x and y, and set
  148. *     up the start points in  X  and  Y.
  149. *
  150.       IF( LSAME( TRANS, 'N' ) )THEN
  151.          LENX = N
  152.          LENY = M
  153.       ELSE
  154.          LENX = M
  155.          LENY = N
  156.       END IF
  157.       IF( INCX.GT.0 )THEN
  158.          KX = 1
  159.       ELSE
  160.          KX = 1 - ( LENX - 1 )*INCX
  161.       END IF
  162.       IF( INCY.GT.0 )THEN
  163.          KY = 1
  164.       ELSE
  165.          KY = 1 - ( LENY - 1 )*INCY
  166.       END IF
  167. *
  168. *     Start the operations. In this version the elements of A are
  169. *     accessed sequentially with one pass through A.
  170. *
  171. *     First form  y := beta*y.
  172. *
  173.       IF( BETA.NE.ONE )THEN
  174.          IF( INCY.EQ.1 )THEN
  175.             IF( BETA.EQ.ZERO )THEN
  176.                DO 10, I = 1, LENY
  177.                   Y( I ) = ZERO
  178.    10          CONTINUE
  179.             ELSE
  180.                DO 20, I = 1, LENY
  181.                   Y( I ) = BETA*Y( I )
  182.    20          CONTINUE
  183.             END IF
  184.          ELSE
  185.             IY = KY
  186.             IF( BETA.EQ.ZERO )THEN
  187.                DO 30, I = 1, LENY
  188.                   Y( IY ) = ZERO
  189.                   IY      = IY   + INCY
  190.    30          CONTINUE
  191.             ELSE
  192.                DO 40, I = 1, LENY
  193.                   Y( IY ) = BETA*Y( IY )
  194.                   IY      = IY           + INCY
  195.    40          CONTINUE
  196.             END IF
  197.          END IF
  198.       END IF
  199.       IF( ALPHA.EQ.ZERO )
  200.      $   RETURN
  201.       IF( LSAME( TRANS, 'N' ) )THEN
  202. *
  203. *        Form  y := alpha*A*x + y.
  204. *
  205.          JX = KX
  206.          IF( INCY.EQ.1 )THEN
  207.             DO 60, J = 1, N
  208.                IF( X( JX ).NE.ZERO )THEN
  209.                   TEMP = ALPHA*X( JX )
  210.                   DO 50, I = 1, M
  211.                      Y( I ) = Y( I ) + TEMP*A( I, J )
  212.    50             CONTINUE
  213.                END IF
  214.                JX = JX + INCX
  215.    60       CONTINUE
  216.          ELSE
  217.             DO 80, J = 1, N
  218.                IF( X( JX ).NE.ZERO )THEN
  219.                   TEMP = ALPHA*X( JX )
  220.                   IY   = KY
  221.                   DO 70, I = 1, M
  222.                      Y( IY ) = Y( IY ) + TEMP*A( I, J )
  223.                      IY      = IY      + INCY
  224.    70             CONTINUE
  225.                END IF
  226.                JX = JX + INCX
  227.    80       CONTINUE
  228.          END IF
  229.       ELSE
  230. *
  231. *        Form  y := alpha*A'*x + y.
  232. *
  233.          JY = KY
  234.          IF( INCX.EQ.1 )THEN
  235.             DO 100, J = 1, N
  236.                TEMP = ZERO
  237.                DO 90, I = 1, M
  238.                   TEMP = TEMP + A( I, J )*X( I )
  239.    90          CONTINUE
  240.                Y( JY ) = Y( JY ) + ALPHA*TEMP
  241.                JY      = JY      + INCY
  242.   100       CONTINUE
  243.          ELSE
  244.             DO 120, J = 1, N
  245.                TEMP = ZERO
  246.                IX   = KX
  247.                DO 110, I = 1, M
  248.                   TEMP = TEMP + A( I, J )*X( IX )
  249.                   IX   = IX   + INCX
  250.   110          CONTINUE
  251.                Y( JY ) = Y( JY ) + ALPHA*TEMP
  252.                JY      = JY      + INCY
  253.   120       CONTINUE
  254.          END IF
  255.       END IF
  256. *
  257.       RETURN
  258. *
  259. *     End of DGEMV .
  260. *
  261.       END
  262.